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Abstract. Computer-assisted prostate biopsies became a very active 
research area during the last years. Prostate tracking makes it possi- 
ble to overcome several drawbacks of the current standard transrectal 
ultrasound (TRUS) biopsy procedure, namely the insufficient targeting 
accuracy which may lead to a biopsy distribution of poor quality, the very 
approximate knowledge about the actual location of the sampled tissues 
which makes it difficult to implement focal therapy strategies based on 
biopsy results, and finally the difficulty to precisely reach non-ultrasound 
(US) targets stemming from different modalities, statistical atlases or 
previous biopsy series. The prostate tracking systems presented so far 
are limited to rigid transformation tracking. However, the gland can get 
considerably deformed during the intervention because of US probe pres- 
sure and patient movements. We propose to use 3D US combined with 
image-based elastic registration to estimate these deformations. A fast 
elastic registration algorithm that copes with the frequently occurring 
US shadows is presented. A patient cohort study was performed, which 
yielded a statistically significant in-vivo accuracy of 0.83±0.54mm. 



1 Introduction 

Prostate biopsies are the only definitive way to confirm a prostate cancer hy- 
pothesis. The current clinical standard is to perform prostate biopsies under 2D 
TRUS control. The US probe is equipped with a needle guide for transrectal 
access of the prostate. The guide aligns the needle trajectory with the US image 
plane, which makes it possible to visualize the trajectory on the image for needle 
placement control. Unfortunately, in particular mid- and early-stage carcinoma 
are mostly isoechogenic, i.e. not visible in US images, which makes it necessary 
to sample the gland according to a systematic pattern. It is common to acquire 
10 to 12 systematically distributed biopsies, the standard pattern taking roughly 
into account that most tumors (70%) develop in the peripheral zone of the gland. 
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The current standard biopsy procedure has several shortcomings: first, it is 
difficult for the clinician to reach systematic targets accurately because he has 
to move the probe continuously to place the needle; a constant visual reference 
is hence lacking. Second, performing a non-exhaustive systematic search for an 
invisible target implies that the target can be missed. Negative results leave the 
clinician in a dilemma when the cancer hypothesis cannot be discarded: his only 
option is to repeat the biopsy series. Furthermore, the location of the acquired 
samples with respect to the patient anatomy is only very approximately known 
after the intervention. Uncertainty about tumor location is the principal reason 
why prostate therapy is in general radical. 

In order to address these issues, Baumann et al. ^ and Xu et al. [2, simul- 
taneously proposed to acquire a US volume before the intervention and to use 
it as anatomical reference. The stream of US control images acquired during the 
intervention is then registered with the reference volume, which allows to project 
targets defined in the reference volume into the control images, and, conversely, 
the biopsy trajectory, known in control image space, into the reference volume. 
This technique makes it possible to improve biopsy distribution accuracy by 
showing the current trajectory in a fixed reference together with the trajectories 
of previously acquired biopsies, to aim targets defined in the reference volume 
during a planning phase, and to know the precise biopsy positions after the 
intervention. Non-US targets could originate from suspicious lesions in MR vol- 
umes that are then multi-modally registered with the US reference volume. It 
is also possible to derive targets from more sophisticated statistical atlases [3] 
or, in the case of repeated biopsies, they could consist of previously unsampled 
regions. After the intervention, the biopsy trajectories in the reference volume 
can be combined with the histological results and used for therapy planning. 

Xu et al. p] acquire a freehand 3D US volume and use 2D control images 
during the intervention. The 2D control images are tracked in operating room 
space with a magnetic sensor on the probe. In a second step, image-based regis- 
tration is performed to compensate for small organ and patient movements. A 
similar approach was proposed by Bax et al. 4 , who use an articulated arm for 
2D US beam tracking. Bax does not compensate for patient and gland move- 
ments. However, pain-related pelvis movements are frequent, since the patient 
is not under total anesthesia. In that case, both methods risk to loose track of 
the gland because the US beam is tracked in operating room space and not in 
organ space, and a new reference volume has to be acquired. Baumann et al. 
address this draw-back by using 3D US to obtain richer control images during 
the intervention . Instead of using a US beam tracking device to initialize local 
image-based registration, they propose a kinematic model of endorectal probe 
movements to compute anatomically plausible positions of the US beam with 
respect to the gland, which is unaffected by patient movements. 

However, probe movements during needle placement continuously deform the 
gland. Deformations are strongest near the probe head and can reach 3 to 6 mm. 
They cannot be estimated with the presented systems. To address this issue, we 
extend the 3D US rigid registration approach presented in [I] by adding a de- 



formation estimation step to the registration pipeline. 3D US control images 
provide the information required to estimate the deformation with acceptable 
precision and accuracy. Inverse consistency and linear elasticity are used as de- 
formation priors. A novel image distance measure capable of dealing with local 
intensity shifts, frequent in US images, is presented. The clinical accuracy of 
the presented algorithm is evaluated on a large number of patient data acquired 
during prostate biopsy sessions. 

2 Method 

Non-linear registration of 3D US image streams is currently the most promising 
approach to perform organ tracking with deformation estimation. The principal 
challenges of image-based tracking systems are robustness and computational ef- 
ficiency. A technique to achieve both goals are coarse-to-fine registration strate- 
gies that successively increase the degrees of freedom (DOF) of the transforma- 
tion space and the image resolution. In this paper, we add an elastic registration 
step to the 3-step coarse-to-fine rigid registration pipeline that we proposed in 
[T]. The resulting pipeline is illustrated in Fig.[TJ 
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Fig. 1. Registration pipeline. The dimensionality of the transformation space 
and the image resolution are successively increased. 



2.1 Ftamework for non-linear registration 

Image-based deformation estimation can be formulated as an optimization pro- 
cess of a local distance measure. Let /i , /2 : — *■ R be images, 93 : ^ the 
deformation function and the functional ^[/i, 1^9] a measure of the distance 
between Ji and I2 o (p. In contrast to parametric approaches that use basis func- 
tions to build the deformation function, we will follow a variational approach and 
define ip{x) = x + u{x), where m : M'^ — > IR^ is assumed to be a diffeomorphism. 
The deformation could then be estimated by solving the optimization problem 



(/3* = arg min (f [/i , /2 ; ) , 

V 



(1) 



where the registration energy £ simply corresponds to T). Straightforward min- 
imization of a distance measure yields, however, in general poor results due to 
countless local minima, in particular in presence of noise, partial object occlusion 
and other imperfections in the image data. Unfortunately, US is a particularly 
noisy modality, which makes 3D US based deformation estimation vulnerable to 
local misregistrations. This problem can be addressed by integration of a pri- 
ori models of the expected deformation. This can be done implicitly by adding 
further energy terms to the objective function. In this work, inverse consistency 
and elastic regularization energies are added. 

2.2 Inverse Consistency Constraints 

In non-linear image registration, the forward estimation that minimizes 5 [/i , /2 ; V 
does in general not yield the inverse of the backward estimation that minimizes 
£[l2,Ii; i.e. ip o ip Id with /c? : R'^ — > M^, x ^ x. Introduction of Zhang's 
inverse consistency constraint ^ 

I[^P;ip] = [ 11^-0^-/^1123 dx (2) 

as additional energy penalizes solutions that lead to inconsistent inverse trans- 
formations, where f2 C K'^ is the registration domain in image space. Estimation 
of the forward and the backward deformations is coupled by an alternating it- 
erative optimization 

^'=+1 = argmin(f[/i,/2;(^] , (3) 

V^'^+i = argmin(£[/2,/i; V] ■ (4) 

Concurrent estimation with mutual correction reduces the risk of local misreg- 
istrations. 

2.3 Elastic Regularization 

The deformation of the prostate caused by probe pressure is fully elastic, which 
justifies the introduction of the linearized elastic potential [6j 



£[ip] =S[u + Id] 




as additional energy, where A and fi are the Lame coefficients. 
2.4 Image Distance Measure 

The image distance measure is the driving energy of the optimization process. 
Experiments on patient data have shown that the sum of squared distances (SSD) 



is a poor distance measure for deformation estimation on noisy US images. Local 
intensity changes are frequent due to changing US beam angles with respect to 
the tissues and probe pressure variations. The more robust Pearson correlation 
coefhcient (CC) requires the evaluation of a large neighborhood of every voxel 
pair to yield statistically significant results, which is incompatible with deep 
multi-resolution approaches that operate on very coarse levels. 

We hence prefer an intermediate correlation model that filters low-frequency 
intensity shifts, i.e. we assume that Ii = I2 o (p + b, where (p is the physical 
solution of the registration problem, and where 6 : M'^ ^ models a local 
intensity shift. The shift is estimated by 

b^[ip]{x)^{h~i2 0ip)*g,{x) (6) 

where ^ : E'^ — > R is a Gaussian with standard deviation a. The image distance 
energy is then 

Vih.h;^] = I [h[x)~h{v{x))~b''M{x)f dx. (7) 
J n 

The standard deviation a controls the frequency range of the high-pass fil- 
ter. If a gets smaller, the cropped frequency range gets larger, and registration 
convergence rate decreases and may even stall if only high frequency noise like 
speckle is left. When used with a multi-resolution solver on a Gaussian pyramid 
(cf. next section) , which implicitly performs a low-pass filtering of the intensity 
variations on coarse resolutions, this approach transforms to a band-pass filter- 
ing on varying frequency bands. In this configuration it is sufficient to chose 
relatively small standard deviations without risking registration inefficiencies. 



2.5 Solver 

Combination of the energy terms yields the alternating system 

^* = argmin {V[h,h. v\ + £M +1^. f]) , (8) 
iP* = a.rgmm{V[l2,Ii,'ip]+£[tp] +I[^;V'])- (9) 

An iterative two-step minimization scheme is used to solve both objective func- 
tions. The Euler-Lagrange equations of Eqn. [5] and [S] are rewritten as a fixed 
point iteration 

/e+l k 

"'^ + /p[/i,/2;/] + /x[V'';/], (10) 



At 

'"^ =c[^^] + h[h,iu^''] + h[v^;^% (11) 



At 

where t G M controls the discretization granularity, and with the elliptic partial 
differential operator 



C[ip\ = C[u + Id] = ^xAu +{X + /z)Vdiv u, 



(12) 



which is obtained from the Gateaux-derivative of £[(p]|6j. The Gateaux deriva- 
tives of the energy term V at tp yields the force term 



/a; = (/i - /2 - * Ga) ■ (V/2) o ^. 



(13) 



and for T we get 



/r[V'; 'f] = (V-fA) o 1^ ■ {iIj o 1^ - Id). 



(14) 



An iterative algorithm is used to estimate the displacement fields: 
1: while not converged do 
2: compute fv[hj2\ V^] and fi[tp^] (fi''] 
3: compute /p [h, h ; ^''] and i^^] 
4: solve Eqn. [To] for 
5: solve Egn. [TTjfor ^^+^ 
6: end while 

The forces are hence considered as constants for the resolution of the PDEs 
[TUl and [TTl and the forward and the backward estimation correct themselves 
mutually at each force update. The PDEs are solved using Red-Black Gauss- 
Seidel relaxation. Convergence is achieved if the difference of the L2-norm of 
the total forces between two iterations is below a threshold for both the forward 
and the backward estimation (oscillatory states are detected). The algorithm is 
executed on various resolution levels of a Gaussian image pyramid [T] using the 
full multigrid strategy [7]- Note that the algorithm derives from the multigrid 
scheme by iterating until convergence on every grid level. This is necessary since 
the relaxation is performed on fractional forces. Fixed edges and bending side 
walls are used as border conditions [5]. The elasticity parameters are chosen 
such that Poisson's coefBcient is zero, hence maximizing compressibility to allow 
compensation of local model inadequacies. Young's modulus is interpreted as a 
free variable in function of Poisson's coefficient and the PDE discretization At 
since it has no physical meaning in image registration. The forces are capped to a 
maximum length which makes it possible to control the maximum contributions 
per iteration to the displacement field via At. Limiting the contributions to 
less than 0.5 voxel side lengths ensures that the algorithm does not 'jump' over 
intensity barriers during optimization. 

3 Experiments 

The framework was validated on 278 registrations of 295 US volumes from 17 
patients. The 17 reference images were acquired shortly before the intervention, 
and the tracking images were acquired after a biopsy shot. The clinical protocol 
was approved by the ethical committee of the XXX hospital. Town, Country, and 
all patients consented to participate to the study. The images were acquired with 
a GE Voluson and a RIC5-9 endorectal US probe. The algorithms were executed 
on a 4-core 2.6Ghz processor. In order to provide a reference gold standard for 
the evaluation of registration accuracy, experts manually segmented 467 point 
fiducials that were clearly identifiable on multiple images (e.g. calcifications and 



cysts). The distances between fiducial pairs were measured after registration 
to estimate the local accuracy. Note that the unavoidable segmentation error 
increases the measured error in average; this approach hence underestimates 
accuracy. Accuracy was computed for all registrations that were qualified as valid 
by experts after visual inspection, which represent 97,8% of the registrations. The 
results for both rigid and elastic registration are given in Tab. [1] and a visual 
illustration of the registration performance is given in Fig. [21 Fig. [3] shows 3D 
biopsy maps created with our biopsy tracking system. 





mean 


standard 


max 


execution 




distance 


deviation 


distance 


time (mean) 


unregistered 


13.76 mm 


7.89 mm 


51.61 mm 




rigid 


1.33 mm 


0.85 mm 


4.19 mm 


2.1 s 


elastic 


0.83 mm 


0.54 mm 


4.14 mm 


6.8 s 



Table 1. Accuracy study. 




Fig. 2. Fig. (a) shows a prostate volume with calcifications [1,2]. Fig. (b) shows a 
second volume after rigid registration; low probe pressure led to the low contrast 
zone [3]. Fig. (c) shows the 3D elastic registration with standard SSD; the whole 
prostate is dragged towards zone [3] . Fig. (d) shows the 3D intensity shift filtered, 
inverse consistent elastic registration; the strong intensity differences between 
both volumes are correctly handled, the calcifications make appearance at the 
correct position (best viewed in PDF with zoom). 



4 Discussion and conclusion 

Deformation estimation achieves an overall accuracy of at least 0.83±0.54 mm on 
real patient data. This corresponds to an error reduction of 40% when compared 
to rigid 3D-3D registration. The average computation time of the registration was 
only 6.8s. We are confident that the algorithm can be accelerated to below Is on 
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Fig. 3. Tissue deformation corrected 3D biopsy maps in reference space. 



the same machine with simple optimization and parallehzation techniques, which 
is sufficient for assisted needle placement. With specialized standard hardware 
(GPUs), at least 5Hz should be feasible. 

Biopsy tracking systems potentially add significant clinical value to prostate 
cancer diagnosis and therapy planning. Immediate advantages are the possibility 
to avoid resampling of already biopsied tissues when repeating a biopsy series, 
interventional quality control of the biopsy distribution (e.g. detection of un- 
sampled areas) and computer-assisted guidance to non-systematic targets. The 
latter could for example be identified on MR/spectroMR images of the gland. 
Moreover, the improved knowledge about the biopsy and thus the cancer posi- 
tion could be used to implement focal therapy strategies for prostate cancer. 3D 
US based elastic tracking can provide the precision required for such therapeutic 
applications. 
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